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1.  Introduction 


K-dimensional  (KD)  trees  use  binary  space-partitioning  algorithms  to  organize  and  store  data 
points  in  K-dimensional  space/  They  are  particularly  useful  for  efficient  nearest-neighbor  search 
algorithms.  This  report  presents  a  set  of  functions,  written  in  C-i-i-,  that  is  designed  to  be  used  to 
create,  search,  and  delete  KD  trees.  All  of  the  functions  are  based  on  recursive  algorithms.  Tests 
for  measuring  function  performance  are  included,  as  are  examples  for  creating  Voronoi 
diagrams. 

The  functions  that  are  described  in  this  report  have  been  grouped  into  the  yKDTree  namespace, 
which  is  summarized  at  the  end  of  this  report.  The  yKDTree  namespace  relies  exclusively  on 
standard  C-i-i-  operations.  However,  example  code  that  is  included  in  this  report  makes  use  of  the 
yRandom  namespace  for  generating  pseudorandom  numbers  and  the  yBmp  namespace  for 
creating  Voronoi  diagrams. 


2.  Sorting  Tables  —  the  ColumnSortO  Function 


The  ColumnSortO  function  uses  a  stable  merge-sort  algorithm  to  sort  tabulated  data  into 
ascending  order.  The  data  must  be  stored  in  a  2-index  array,  where  the  indices  are  arbitrarily 
referred  to  here  as  rows  (first  index)  and  columns  (second  index). 

The  ColumnSortO  function  is  included  in  the  yKDTree  namespace  as  a  helper  function  for  the 
NewTreeO  function,  which  is  described  in  section  4.  However,  the  ColumnSortO  function  can 
also  be  useful  on  its  own  when  an  efficient  means  of  sorting  tables  by  columns  is  desired. 

2.1  ColumnSortO  Code 

template<class  T>voici  ColumnSort(//<=SORT  A  TABLE  BY  COLUMNS  (ASCENDING  ORDER) 


T**a,T**b,//< - POINTERS  TO  STARTING  &  ENDING  ROWS  (FIRST  INDEX) 

T**t,//< - TEMPORARY  STORAGE  (SIZE  >=  b-a) 

unsigned  c){//< - THE  COLUMN  (2ND  INDEX)  TO  SORT  ON 

if (b-a<2)return; 


T**m=a+(b-a)/2;/*->*/ColumnSort(a,m,t,c),ColumnSort(m,b,t,c); 

for(T**u=m, **l=a, ** j=t; j<b-a+t; )* j++=u<b&&(l>=m | |(*l)[c]>(*u)[c] ) ?*u++: *1++; 

while (a<b)*a++=*t++; 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~15DUL2014 - 
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2.2  ColumnSortO  Parameters 


a  a  points  to  the  first  row  that  will  be  included  in  the  sort. 

b  b  points  to  one  past  the  last  row  that  will  be  included  in  the  sort. 

t  t  points  to  temporary  storage  for  the  ColumnSortO  function,  t  must  point  to  a 

array  that  is  capable  of  storing  at  least  b-a  elements. 

c  c  specifies  the  column  on  which  the  sort  will  be  based. 

2.3  ColumnSortO  Example 

The  following  example  uses  the  ColumnSortO  function  to  sort  a  2-column  table,  first  by  the 
second  column,  then  by  the  first. 

#include  <cstdio>// . printf() 

#include  "y_kd_tree . h"// . yKDT ree 

int  main(){//<====================A  SIMPLE  EXAMPLE  FOR  THE  ColumnSortO  FUNCTION 

int  n[9],A[9][2]={0,5,0,7,l,9,0,8,l,6,0,7,l,4,0,0,0,3},*B[9]; 
for(int  i=0;i<9;++i)B[i]=A[i]; 

yKDTree: :ColumnSort(B,B+9,T,l), yKDT ree: :ColumnSort(B,B+9,T,0); 

printfC  UNSORTED  |  SORTED\n - | - \n"); 

for(int  i=0;i<9;++i) 

printfC  %d  ,  %d  |  %d  ,  %d\n"CA[i]  ,A[i]  [1]  CB[i] ,  B[i]  [1] ) ; 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  D^l  S  D  U  L  2  0 


OUTPUT: 


UNSORTED 

!  1 

SORTED 

0 

> 

5 

1 

1 

0  ,  0 

0 

) 

7 

1 

0  ,  3 

1 

> 

9 

1 

0  ,  5 

0 

) 

8 

1 

0  ,  7 

1 

) 

6 

1 

0  ,  7 

0 

> 

7 

1 

0  ,  8 

1 

) 

4 

1 

1  ,  4 

0 

> 

0 

1 

1  ,  6 

0 

i 

3 

1 

1  ,  9 

2.4  ColumnSortO  Performance 

The  following  code  measures  the  performance  of  the  ColumnSortO  function  and  compares  it 
with  the  performance  of  the  stable_sortO  function.  The  yRandom  namespace  is  used  to  generate 
pseudorandom  numbers  for  the  test.  Time  measurements  represent  the  total  time  required  to 
perform  10®  sorts  for  tables  with  2"  rows,  where  n  varies  from  1  to  14. 
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The  output  was  generated  by  compiling  the  code  using  Microsoft’s  Visual  Studio  C++  2010 
Express  compiler,  with  the  output  set  to  “release”  mode.  For  this  scenario,  the  ColumnSort() 
function  outperforms  the  built-in  stable_sort()  function  (Fig.  1). 


#include 

#include 


#include  <algorithm>// . stable_sort( ) 

#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) ,CLOCKS_PER_SEC 

y_kd_tree .  h"// . yKDT  ree 

y_random.h"// . yRandom 

inline  bool  Compare(//<======HELPER  FUNCTION  FOR  THE  std : : stable_sort( )  FUNCTION 

double*a,double*b){//< - POINTERS  TO  COMPARISON  ROWS 

return  a[0]<b[0];// . note  column  number  is  fixed  at  0 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T^U  PDAT  EO'^ISDUL^G 

int  main(){//<==============MEASURE  THE  PERFORMANCE  OF  THE  ColumnSort()  FUNCTION 

const  int  N=1<<14,M=1000000;//. . . .max  #  of  rows,  number  of  iterations  per  test 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 

double  s,t,*T[N],*A[N];/*<-*/for(int  i=0;i<N;++i)A[i]=new  double[l]; 
printf("  I  std ; ; stable_sort( )  |  ColumnSort( ) \n" 

"  row  I - I - \n" 

"  count  I  time  |  Z[m/2][0]  |  time  |  Z[m/2][0]\n" 

"  I  (s)  I  avg.  I  (s)  I  avg.\n" 

"  - I - I - I - I - \n");//table  header 


for(int  m=2;m<=N;m*=2){ 

s=0,t=clock(), yRandom: : Initialize(1, 1) ; // . begin  stable_sort( )  test 

for(int  k=0; k<M;++k){ 

for (int  i=0; i<m;++i)A[i] [0]=yRandom: :RandU(I,0,l); 
std : : stable_sort ( A, A+m, Compare ) , s+=A[m/2] [0] ; } 
printf("%7d  |%9.3f  |%10.7f  | ",m, (clock()-t)/CLOCKS_PER_SEC,s/M); 

s=0,t=clock( ), yRandom: : Initialize(1, 1); // . begin  ColumnSort( )  test 

for(int  k=0; k<M;++k){ 

for (int  i=0; i<m;++i)A[i] [0]=yRandom: :RandU(I,0,l); 
yKDTree: :ColumnSort(A,A+m,T,0) , s+=A[m/2] [0] ; } 
printf("%9.3f  | %10.7f \n" , (clock( ) -t)/CLOCKS_PER_SEC, s/M) ; } 
for (int  i=0;i<N;++i)delete[ ]A[i] ; 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  '**'U  PDAT  ED'^TSDUL20 
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OUTPUT: 


row 

count 

1  std : : stable  sort( )  | 

1  1 

ColumnSortO 

1 

1  time 

1  (s) 

Z[m/2][0]  i 
avg.  1 

time 

(s) 

Z[m/2][0] 

avg. 

2 

1  0.046 

0.6667224  | 

0.032 

0.6667224 

4 

1  0.078 

0.6002211  1 

0.093 

0.6002211 

8 

1  0.219 

0.5556498  j 

0.249 

0.5556498 

16 

1  0.546 

0.5293284  j 

0.577 

0.5293284 

32 

1  1.316 

0.5150998  1 

1.263 

0.5150998 

64 

1  3.401 

0.5077139  1 

2.886 

0.5077139 

128 

1  7.675 

0.5038248  1 

6.646 

0.5038248 

256 

1  17.456 

0.5019407  1 

14.682 

0.5019407 

512 

1  39.059 

0.5009805  1 

31.961 

0.5009805 

1024 

1  83.569 

0.5004688  1 

68.828 

0.5004688 

2048 

1  186.202 

0.5002375  1 

150.859 

0.5002375 

4096 

1  409.029 

0.5001205  1 

324.842 

0.5001205 

8192 

1  890.498 

0.5000599  1 

697.758 

0.5000599 

16384 

2061.582 

0.5000335 

1537.217 

0.5000335 

2500 


0  2048  4096  6144  8192  10240  12288  14336  16384 


Row  Count 


Fig.  1  ColumnSortO  performance  compared  with  stable_sort()  performance 


3.  KD-Tree  Nodes  —  the  NODE  Struct 


NODE  structs  can  be  used  to  store  the  nodes  that  make  up  KD  trees.  Typieally,  NODE  struets 
are  ereated  using  the  NewTree()  funetion  (see  seetion  4)  and  deleted  using  the  DeleteTree() 
funetion  (see  seetion  5). 
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3.1  NODE  Code 


template<class  T>struct  lsiODE{//<================== 

T*r;//< - 

NODE*a,*b;//< - 

int  k;//< - 


==============A  KD-TREE  NODE 

-POINTER  TO  A  ROW  IN  A  TABLE 
- POINTERS  TO  SUBNODES 

- node  layer 

LAST~UPDATED~15DUL2014 - 


3.2  NODE  Parameters 

r  r  points  to  a  row  in  a  table, 

a  a  points  to  a  subnode, 

b  b  points  to  a  subnode, 

k  k  is  used  to  identify  a  node’s  layer. 


4.  Creating  KD  Trees  —  the  NewTreeO  Function 


If  a  2-index  array  is  used  to  store  sortable  tabulated  data  in  the  format  that  is  shown  in  Fig.  2, 
then  the  NewTree()  function  can  be  used  to  create  a  KD  tree  for  the  tabulated  data.  In  Fig.  2, 
values  for  both  indepenent  and  dependent  variables  are  stored  in  array  A  .  Independent  variables 
are  stored  in  columns  with  subscripted-  X  headers,  while  dependent  variables  are  stored  in 
columns  with  subscripted- F  headers.  Each  row  represents  a  single  data  point. 


Xo 

■ 

••  ■ 

A 

A 

0 

A),o 

A,i 

A).k 

A.x-i 

A.x 

A,x+i 

1 

A,o 

A,i  ■ 

••  At  ■ 

A./r-i 

A.x 

A 

^l,X+l 

i 

A,o 

At  ■ 

••  At  ■ 

Ax-i 

Ax 

A,x+i 

m-l 

A 

^m-1,0 

A 

■  ■A 

A 

A 

^m-\.K+\ 

Fig.  2  Tabulated  data  stored  in  a  2-index  array 
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Note  that  the  NewTree()  function  uses  the  “new”  command  to  allocate  memory  for  nodes.  To 
avoid  memory  leaks,  the  DeleteTree()  function  (see  section  5)  should  be  used  to  deallocate 
memory  when  a  KD  tree  is  no  longer  needed. 

4.1  NewTreeO  Code 

template<class  T>NODE<T>*NewTree(//<==========================CREATE  A  KD  TREE 

T**a,T**b,//< - POINTERS  TO  STARTING  AND  ENDING  ROWS 

int  K,//< - NUMBER  OF  INDEPENDENT  DIMENSIONS 

int  k=0){//< - NODE  LAYER  (SET  BY  RECURSIVE  ALGORITHM) 

T**t=new  T* [b-a] ; /*->*/ColumnSort(a, b,t, k) ; /*&*/cielete[ ]t; 

T**m=a+(b-a-l)/2; 

NODE<T>*N=new  NODE<T>;/*<-*/N->r=*m,N->a=a-m?NewTree(a,m,K, (k+l)%K) :0, 
N->b=b-m-l?NewTree(m+l,b,K, (k+l)%K) :0,N->k=k; 
return  N; 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~15DUL2014 - 


4.2  NewTreeO  Parameters 

a  a  points  to  the  first  row  of  the  table  that  will  be  included  in  the  new  KD  tree. 

b  b  points  to  one  past  the  last  row  of  the  table  that  will  be  included  in  the  new  KD 

tree. 

K  K  specifies  the  number  of  independent  dimensions  that  are  included  in  the  table 

that  is  specified  by  a  and  b. 

k  k  specifies  the  layer  of  the  current  node  that  is  being  created  by  the  NewTreeO 

function,  k  is  set  automatically,  either  by  the  default  value,  or  by  the  NewTreeO 
function  when  it  calls  itself  recursively. 

4.3  NewTreeO  Return  Value 

The  NewTreeO  function  returns  a  pointer  to  the  root  node  of  a  newly  created  KD  tree. 

4.4  NewTreeO  Example  —  Creating  a  Simple  KD  Tree 

The  following  example  code  uses  the  NewTreeO  function  to  create  a  simple  KD  tree,  with 

K  =  2  independent  dimensions,  then  prints  the  nodes.  Nodes  represented  by  (X,X)  are  empty. 
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#include  <cstdio>// . printf() 

#include  "y_kd_tree . h"// . yKDT ree 

inline  void  PrintNode(//<=======================PRINTS  THE  COORDINATES  OF  A  NODE 

yKDTree:  :NODE<int>*N){//< - A  NODE 

N?printf("(%d,%d)  ",N->r[0],N->r[l]) :printf("(X,X)  "); 

^11  E  N  AUT  ^)GI^AI  L  •  LAS  T  PDAT  E  D^l  S  D  U  L  2  0 

int  main(){//<=====================A  SIMPLE  EXAMPLE  FOR  THE  NewKDTree()  FUNCTION 

int*A[5],B[5] [3]={4,1,0  ,  4,3,1  ,  6,2,2  ,  2,4,3  ,  8,4,4}; 
printf( "TABULATED  DATA:\n"); 

for (int  i=0; i<5;++i)printf ( "%13d,%d,%d\n'', B[i] [0],B[i][l],B[i][2]); 

printf("\n\nKD  TREE:\n"); 

for(int  i=0;i<5;++i)A[i]=B[i]; 

yKDTree: :NODE<int>*N=yKDTree: :NewTree(A,A+5,2); 

printf("  " ) , PrintNode(N) ; // . root  node  (k=0) 

printf^'Xn  /  \\\n  "); 

PrintNode(N->a),printf("  ");PrintNode(N->b);//. . . .2nd-level  nodes  (k=l) 

printf("\n  /  \\  /  \\\n"); 

PrintNode(N->a->a),PrintNode(N->a->b),PrintNode(N->b->a),PrintNode(N->b->b); 
printf("\n\n"), yKDT ree: :DeleteTree(N); 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  '**'U  PDAT  ED'^TSDUL20 


OUTPUT: 


TABULATED  DATA: 

4,1,0 

4,3,1 

6,2,2 

2,4,3 

8,4,4 

KD  TREE: 

(4,3) 

/  \ 

(4,1) 

(6,2) 

/  \ 

/  \ 

(X,X)  (2,4)  (X,X)  (8,4) 

5.  Deleting  KD  Trees  —  the  DeleteTreeO  Function 

The  DeleteTreeO  function  can  be  used  to  delete  a  KD  tree  that  has  been  created  using  the 
NewTreeO  function. 

5.1  DeleteTreeO  Code 
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template<class  T>void  DeleteTree(//<==========================DELETE  A  KD  TREE 

NODE<T>*N){//< - THE  ROOT  NODE  OF  A  KD  TREE 

if ( !N) return; /*else*/DeleteTree(N->a),DeleteTree(N->b), delete  N; 

}// -  ~YAGENAUT@GMAI L . 


5.2  DeleteTreeO  Parameters 

N  N  points  to  the  root  node  of  a  KD  tree. 

6.  Searching  KD  Trees  —  the  NNSearch()  Function 


The  NNSearchO  function  can  be  used  to  search  a  KD  tree  for  ,  a  row  in  a  table  that  specifies  a 
location  that  minimizes  the  distance  to  X  ,  a  user  specified  location  where 

(1) 

Thus,  A. is  defined  to  be  a  row  for  which  S  in  Eq.  2  is  minimized. 

S  =  T{x,-aJ\  (2) 

k=0 


Note  that  A,  may  or  may  not  be  unique. 

6.1  NNSearchO  Code 

template<class  T, class  U>void  NNSearch(//<=============NEAREST-NEIGHBOR  SEARCH 

const  NODE<T>*N,//< - THE  ROOT  NODE  OF  A  KD  TREE 


const  T*X,//< - POINTER  TO  SEARCH  COORDINATES 

T*&r,//< - POINTER  TO  NEAREST-NEIGHBOR  ROW  (CALCULATED) 

U&S,//< - SQUARED  DISTANCE  BETWEEN  COORDINATES  AT  X  AND  r  (CALCULATED) 

int  K){//< - NUMBER  OF  INDEPENDENT  DIMENSIONS 


NODE<T>*a, *b; /*<-*/X[N->k]<N->r[N->k] ?a=N->a,b=N->b: (a=N->b,b=N->a); 
if (a)NNSearch(a,X,  r,S,  K) ; 

U  s=0; /*<-*/for(int  k=0; k<K;++k)s+=(X[k] -N->r[k])*(X[k]-N->r[k]); 
if (s<S)r=N->r,S=s; 

if((s=X[N->k]-N->r[N->k])*s<S&&b)NNSearch(b,X,r,S,K); 

}//  YAGENAUT@GMAI  L  . 


6.2  NNSearchO  Parameters 

N  N  points  to  the  root  node  of  a  KD  tree. 

X  X  points  to  X  ,  the  coordinates  of  the  location  for  the  search.  The  coordinates  that 

are  pointed  to  by  X  must  consist  of  K  values  that  correspond  to  the  K  independent 
variables  that  are  associated  with  the  KD  tree  specified  by  N. 
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r  r  points  to  A. ,  a  row  that  specifies  the  coordinates  of  a  location  for  which  the 

distance  to  the  location  specified  by  X  is  minimized,  r  is  calculated  by  the 
NNSearchO  function. 

S  S  specifies  S  ,  the  squared  distance  between  the  locations  specified  by  X  and  r. 

Although  S  is  calculated  by  the  NNSearch()  function,  S  should  be  initialized  to 
some  value  that  is  larger  than  the  expected  final  value  of  S  (typically  some  very 
large  value). 

K  K  specifies  K ,  the  number  of  independent  dimensions  that  are  included  in  the  KD 

tree. 


6.3  NNSearchO  Example  —  Creating  a  Simple  Voronoi  Diagram 


The  following  example  code  uses  functions  from  the  yBmp  namespace,  along  with  the 
NNSearchO  function  to  create  the  Voronoi  diagram  that  is  presented  in  Fig.  3.  The  black  dots  in 
Fig.  3  show  the  locations  of  the  points  that  were  used  to  create  the  Voronoi  diagram.  The  colored 
sections  represent  sets  of  points  that  have  a  common  nearest  neighbor  among  the  points  that 
make  up  the  KD  tree  (i.e.,  the  black  dots).  For  this  particular  Voronoi  diagram,  the  colors 
themselves  represent  the  row  number  of  the  table  that  was  used  to  create  the  KD  tree  (see  Fig.  4). 


#include  <cmath>// . fabs() 

#include  "y_bmp. h"// . yBmp, <cstring>{memcpy( )} 

#include  "y_kd_tree . h"// . yKDT ree 

inline  void  Rainbow(//<========================================RAINBOW  COLOR  MAP 

unsigned  char  C[3],//< - OUTPUT  COLOR  (CALCULATED) 

double  x,//< - VALUE  FOR  WHICH  A  COLOR  WILL  BE  CALCULATED 

double  min, double  max){//< - MINIMUM  AND  MAXIMUM  SCALED  VALUES 

if (x<min){C[0]=C[l]=C[2]=0; /*&*/return; }// . set  too  small  values  to  black 

if (x>max){C[0]=C[l]=C[2]=255; /*&*/return; }// . set  too  large  values  to  white 

x=(l-(x-min)/(max-min))*8;// . remap  x  to  a  range  of  8  to  0 

C[0]=int((3<x&&x<5| |x>7  ?-fabs(x/2-3)+1.5:5<=x&&x<=7?l:0)*255);// . blue 

C[l]=int((l<x&&x<3i i5<x&&x<7?-fabs(x/2-2)+1.5:3<=x&&x<=5?l:0)*255);// _ green 

C[2]=int((  x<li i3<x&&x<5?-fabs(x/2-l)+1.5:l<=x&&x<=3?l:0)*255);// . red 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  '**'U  PDAT  ED'^ISDUL 

int  main(){//<=========================CREATE  A  VORONOI  DIAGRAM  USING  NNSearchO 

double*A[5],B[5][3]={4,l,0  ,  4,3,1  ,  6,2,2  ,  2,4,3  ,  8,4,4}; 
for(int  i=0;i<5;++i)A[i]=B[i]; 

yKDTree: :NODE<double>*N=yKDTree: :NewTree(A,A+5,2); 

int  n=1000;// . image  size  will  be  2n  x  n  pixels 

unsigned  char*I=yBmp: :NewImage(2*n,n,255),BLACK[3]={0}; 
for(int  p=0;p<2*n;++p)for(int  q=0;q<n;++q){ 
double  X[2]={p*10./(2*n-l),q*5./(n-l)}; 
double  S,*r; /*<-*/yKDTree: :NNSearch(N,X,r,S=lE9,2); 
if (S< .002)memcpy(yBmp: :GetPixel(I, p,q) , BLACK, 3) ; 
else  Rainbow(yBmp: :GetPixel(I, p,q) , r [2] , -1 . 5, 5. 5) ; } 
yBmp : :WriteBmpFile( "voronoi . bmp" , I ) ; 
yKDT  ree : : DeleteT  ree ( N ) , delete [ ] I ; 


}// 


YAGENAUT@GMAIL.COM 


LAST~UPDATED~15DUL2014 
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Fig.  4  Scale  for  Voronoi  diagrams 

6.4  NNSearchO  Performance 

The  following  code  measures  the  performance  of  the  NNSearch()  function  and  compares  it  with 
the  performance  of  the  NNSearchExhaustive()  function,  which  is  defined  in  the  example.  The 
NNSearchExhaustiveO  function  takes  a  brute-force  approach  to  determining  A-  and  S  . 

The  yRandom  namespace  is  used  to  generate  pseudorandom  numbers  for  the  test.  Time 
measurements  represent  the  total  time  required  to  perform  10^  searches  on  tables  with  2”  rows, 
where  n  varies  from  1  to  14. 

When  K  =  2,  the  test  shows  that  for  tables  with  very  few  rows  (somewhere  around  32  or  fewer) 
the  brute-force  method  outperforms  the  KD-tree  method.  Eigure  5  shows  that  as  the  value  of  K 
increases,  the  minimum  number  of  rows  required  for  the  KD-tree  method  to  be  advantageous 
increases  as  well. 
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#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) ,CLOCKS_PER_SEC 

#include  "y_kd_tree . h"// . yKDT ree 

#include  "y_r'andom. h"// . yRandom 

template<class  T>T  NISISearchExhaustive(//<=====EXHAUSTIVE  NEAREST-NEIGHBOR  SEARCH 

T**a,T**b,//< - POINTERS  TO  STARTING  AND  ENDING  ROWS 

T*X,//< - POINTER  TO  SEARCH  COORDINATES 

T*&r,//< - POINTER  TO  NEAREST-NEIGHBOR  ROW  (CALCULATED) 

int  K){//< - NUMBER  OF  INDEPENDENT  DIMENSIONS 

double  S=0; /*<-*/for(int  i=0; i<K;++i)S+=(a[0] [i]-X[i])*(a[0][i]-X[i]); 
for(r=*a;a<b;++a){ 

double  s=0; /*<-*/for(int  i=0; i<K;++i)s+=( (*a) [i]-X[i])*((*a)[i]-X[i]); 
if (s<S)r=*a,S=s; } 
return  S; 

^11  ^'**'^'**'^'**' YAG  E  N  AUT  ^^GI^AI  L  •  LAS  T  '**'U  PDAT  EO'^1SDUL20 

int  main(){//<================MEASURE  THE  PERFORMANCE  OF  THE  NNSearch()  FUNCTION 

const  int  N=1<<14,M=10000000, K=2; 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 
double*A[N] , B[N] [K] ; /*<-*/for(int  i=0; i<N;++i)A[i]=B[i] ; 
for (int  i=0;i<N;++i)for(int  k=0; k<K;++k)A[i] [k]=yRandom: :RandU(I,0,l); 
printf("  I  NNSearchExhaustive( ) I  NN2DInterp( ) \n" 

"  pow  I - I - \n" 

"  count  I  time  |  x  |  time  |  x\n" 

"  I  (s)  I  avg.  I  (s)  I  avg.\n" 

"  - I - I - I - I - \n");//.  .table  header 


for(int  m=2;m<=N;m*=2){ 

double  s=0,t=clock();/*&*/yRandom: :Initialize(I,l);//. . .NNSearchExhaustive( ) 
for(int  k=0; k<M;++k){ 

double  X[K];/*<-*/for(int  i=0;i<K;++i)X[i]=yRandom: : RandU(I,0, 1) ; 
double*r; /*<-*/NNSearchExhaustive(A,A+m,X, r, K) ; 
s+=*r; } 

printf("%7d  |%8.3f  |%9.6f  | ",m, (clock()-t)/CLOCKS_PER_SEC,s/M); 

s=0,t=clock( ), yRandom: : Initialize(1, 1); // . begin  NNSearch()  test 

yKDTree: :NODE<double>*R=yKDTree: :NewTree(A,A+m,K); 
for(int  k=0; k<M;++k){ 

double  X[K];/*<-*/for(int  i=0;i<K;++i)X[i]=yRandom: : RandU(I,0, 1) ; 
double  S=1E9; 

double*r;/*<-*/yKDTree: :NNSearch(R,X, r,S,K); 
s+=*r; } 

yKDTree: :DeleteTree(R); 

printf("%8.3f  | %9.6f \n", (clock( ) -t)/CLOCKS_PER_SEC, s/M) ; } 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  '**'U  PDAT  E  D^l  S  D  U  L  2  0 
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OUTPUT: 


row 

count 


NNSearchExhaustive( ) 


time 


1  (s) 

1 

.  1 

avg.  1 

(s) 

1 

.  1 

avg. 

2 

1  0.320 

1 

1 

0.577847  1 

0.390 

1 

1 

0.577847 

4 

1  0.380 

1 

0.372263  1 

0.646 

1 

0.372263 

8 

1  0.552 

1 

0.354816  1 

0.801 

1 

0.354816 

16 

1  0.710 

1 

0.481116  1 

0.920 

1 

0.481116 

32 

1  1.050 

1 

0.491685  1 

1.200 

1 

0.491685 

64 

1  1.760 

1 

0.499086  1 

1.506 

1 

0.499086 

128 

1  3.104 

1 

0.499204  1 

1.760 

1 

0.499204 

256 

1  5.630 

1 

0.498995  1 

2.060 

1 

0.498995 

512 

1  11.131 

1 

0.499451  1 

2.360 

1 

0.499451 

1024 

1  21.392 

1 

0.499869  1 

2.610 

1 

0.499869 

2048 

1  43.283 

1 

0.499950  1 

2.980 

1 

0.499950 

4096 

1  88.706 

1 

0.499952  1 

3.091 

1 

0.499952 

8192 

1  177.382 

1 

0.499965  1 

3.450 

1 

0.499965 

16384 

382.268 

1 

0.499971 

3.880 

1 

0.499971 

NN2DInterp() 


time 


Fig.  5  Brute-force  and  KD-tree  methods  compared  for  K=2,A,  and  6 
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7.  Example  —  Fermat’s  Spiral 


According  to  Vogel'^  Eqs.  3  and  4,  which  define  a  set  of  points  in  polar  coordinates  that  all  lie  on 
Fermat's  spiral,  can  be  used  to  model  the  patterns  of  seeds  in  sunflowers. 


(3) 


6*  =  (3  -  J5)i7r , 


(4) 


where  R  is  the  radius  of  a  circle  that  contains  all  of  the  points,  m  is  the  total  number  of  points, 
andO<  i  <m. 


The  following  code  uses  functions  from  the  yBmp  namespace,  as  well  as  the  Rainbow()  function 
from  section  6.3,  to  create  three  Voronoi  diagrams  (Fig.  6)  that  are  based  on  Eqs.  3  and  4  (one 
with  m  =100,  one  with  m  =200,  and  one  with  m  =600).  For  each,  R  has  been  chosen  to  be  1/2  the 
width  of  the  image. 


#include  <cmath>// . sqrt() 

#include  "y_bmp. h"// . yBmp, <cstring>{memcpy( )} 

#include  "y_kd_tree . h"// . yKDT ree 

int  main(){//<==================CREATE  VORONOI  DIAGRAMS  BASED  ON  FERMAT'S  SPIRAL 

for(int  m=100,j=0,n=1000; j<3;++j,m*=j+l){//. .. .image  size  will  be  n  x  n  pixels 

double**A=new  double*[m];/*<-*/for(int  i=0;i<m;++i)A[i]=new  double[3]; 
for(int  i=0;i<m;++i){ 

double  R=. 5*sqrt(i/ (m-1. )), theta=(3-sqrt (5. )) *1*3.141592653589793; 
double  x=R*cos ( theta ) ,y=R* sin (theta ) ; 

A[i][0]=x+.5,A[i][l]=y+.5,A[i][2]=i;} 
yKDTree: :NODE<double>*N=yKDTree: : NewTree(A,A+m, 2) ; 
unsigned  char*B=yBmp: :NewImage(n,n,255),BLACK[3]={0}; 
for(int  p=0;p<n;++p)for(int  q=0;q<n;++q){ 
double  X[2]={p/(n-l. ),q/(n-l. )}; 

double  S,*r;/*<-*/yKDTree: :NNSearch(N,X,r,S=lE9,2); 
if (S< . 00002 )memcpy (yBmp: :GetPixel(B, p,q) , BLACK,  3) ; 
else  Rainbow(yBmp: :GetPixel(B, p,q),r[2],0,m);} 
yBmp: :WriteBmpFile( j==0?"spirall . bmp" : j==l?"spiral2 . bmp" : "spiral3.bmp" , B) ; 
for (int  i=0; i<m; ++i) delete [ ]A[i] ; /*&*/yKDTree : :DeleteTree(N) , delete [ ]B; } 

"^11  E  N  AUT  ^)GI^AI  L  •  LAS  T^U  PD  AT  E  D^l  S  D  U  L  2  0 
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Fig.  6  Voronoi  diagrams  of  Fermat’s  Spiral  {m  =  100  left,  m  =  200  center,  m  =  600  right) 


8.  Example  —  Comparing  Average  Distances  to  Sources 


Suppose  that  the  effectiveness  of  some  physical  phenomenon  of  interest  is  reduced  as  the 
distance  from  a  source  increases  (such  as  signal  strength  from  radio  antennas).  If  D(x,  y)  is 

defined  to  be  the  distance  to  the  nearest  source,  then  Eq.  5  can  be  used  to  calculate  the  average 
distance  to  the  nearest  source  ( D  )  for  some  area  of  interest: 

y)dxdy 


where  A  is  the  total  area  over  which  D{x,  y)  is  integrated.  By  comparing  D  values  for  different 
sets  of  points,  the  relative  coverage  effectiveness  between  sets  can  be  evaluated. 

The  following  example  code  uses  functions  from  the  yBmp  namespace,  as  well  as  the  Rainbow() 
function  from  section  6.3,  to  create  3  Voronoi  diagrams  (Fig.  7).  The  first  2  are  based  on  tables 
containing  randomly  selected  points,  while  the  third  is  based  on  a  set  of  points  that  was 
purposely  chosen  to  result  in  a  small  value  for  D  . 
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#include  "y_kcl_tree . h"// . yKDT ree 

#include  "y_bmp. h"// . yBmp, <cstring>{memcpy( )} 

#include  "y_random. h"// . yRandom 

int  main(){//<===========================================COMPARE  AREAL  COVERAGES 

const  int  m=9,  n=1000; // . image  size  will  be  n  x  n  pixels 

double*A[m];/*<-*/for(int  i=0;i<m;++i)A[i]=new  double[3]; 
double  d=l./6; 

double  0[m][2]={3*d,3*d,5*d,3*d,5*d,5*d,3*d,5*d,d,5*d,d,3*d,d,d,3*d,d,5*d,d}; 
unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 
unsigned  char*B=yBmp: :NewImage(n,n,255),BLACK[3]={0}; 
for(int  3=0;3<3;++3){ 
for(int  i=0;i<m;++i) 

A[i] [0]=3==2?O[i] [0] : yRandom: : RandU(I,0, 1) , 

A[i] [l]=3==2?0[i] [1] : yRandom: : RandU(I,0, 1) , 

A[i][2]=i; 

yKDTree: :NODE<double>*N=yKDTree: :NewTree(A,A+m,2); 
double  s=0; 

for(int  p=0;p<n;++p)for(int  q=0;q<n;++q){ 
double  X[2]={p*l./(n-l),q*l./(n-l)}; 
double  S,*r;/*<-*/yKDTree: :NNSearch(N,X,r,S=lE9,2); 
s+=sqrt(S); 

if (S< . 00002 )memcpy(yBmp: :GetPixel(B, p,q) , BLACK, 3) ; 
else  Rainbow(yBmp: :GetPixel(B,p,q),r[2], -l,m);} 
yBmp: :WriteBmpFile( ! 3  ?"coveragel . bmp" : 3==l?"coverage2. bmp" : 

"coverage3 . bmp" , B) ; 

printf("CASE  %d:  D_bar=%8.5f\n",3+l,s/n/n); 
yKDTree: :DeleteTree(N); } 
for (int  i=0;i<m;++i)delete[ ]A[i] ;/*&*/delete[ ]B; 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  D^l  S  D  U  L  2  0 


OUTPUT: 


CASE  1:  D_bar=  0.27243 
CASE  2:  D_bar=  0.17623 
CASE  3:  D  bar=  0.12766 


Fig.  7  Voronoi  diagrams  (case  1  left,  case  2  center,  case  3  right) 
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9.  Example  —  Optimizing  Areal  Coverage 


The  following  example  code  begins  by  recreating  the  “CASE  1”  KD  tree  from  the  example 
presented  in  section  8.  Next,  a  random- walk  method  is  used  to  determine  the  optimal  placement 
of  an  additional  point.  Figure  8  shows  the  original  Voronoi  diagram  compared  with  the  Voronoi 
diagram  with  the  additional  point. 

The  purpose  of  this  example  is  to  show  a  type  of  problem  that  might  benefit  from  the  use  of  KD 
trees.  For  this  particular  case,  since  the  number  of  points  in  the  table  being  searched  is  so  small, 
it  would  likely  have  been  slightly  faster  to  use  a  brute-force  method. 


#include  "y_kd_tree . h"// . yKDT ree 

#include  "y_bmp. h"// . yBmp, <cstring>{memcpy( )} 

#include  "y_random. h"// . yRandom 

int  main(){//<===========================================OPTIMIZE  AREAL  COVERAGE 

double*! [10] , *A[10] ;/*<-*/for(int  i=0;i<10;++i)A[i]=new  double [3] ; 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 

int  n=1000;// . image  size  will  be  n  x  n  pixels 

unsigned  char*B=yBmp: :NewImage(n,n,255),BLACK[3]={0}; 
for(int  i=0;i<10;++i) 

A[i] [0]=yRandom: : RandU(I,0, 1) ,A[i] [l]=yRandom: :RandU(I,0,l),A[i] [2]=i; 
double  s=lE99,x=0,y=0, st,xt,yt; 
for(int  3=0;3<100;++3){ 

for(int  i=0;i<10;++i)T[i]=A[i]; 

st=0,xt=*T[9]=x+yRandom: :RandN(I,0, .l),yt=T[9] [l]=y+yRandom: :RandN(I,0, .1); 
if(xt<0| |xt>l| |yt<0| I yt>l)continue; 
yKDTree:  :NODE<double>*ISI=yKDTree:  :NewTree(T,T+10,2); 
for(int  p=0;p<n;++p)for(int  q=0;q<n;++q){ 
double  X[2]={p*l./(n-l),q*l./(n-l)},D,*r; 
yKDTree:  :NISISearch(N,X,  r,D=lE9, 2) ,  st+=sqrt(D) ;  } 
if (st<s)printf ("3=%2d  ,  D_bar=%f  ,  x=%9.6f  ,  y=%9.6f\n",3, 
(s=st)/n/n,x=xt,y=yt); 
yKDTree: :DeleteTree(N); } 

A[9][0]=x,A[9][l]=y; 

yKDTree: :NODE<double>*N=yKDTree: :NewTree(A,A+10,2); 
for(int  p=0;p<n;++p)for(int  q=0;q<n;++q){ 
double  X[2]={p*l./(n-l),q*l./(n-l)}; 
double  D, *r; /*<-*/yKDTree : :NNSearch(N,X, r,D=lE9,2) ; 
if (D< . 00002 )memcpy(yBmp: :GetPixel(B,p,q),BLACK,3); 
else  Rainbow(yBmp: :GetPixel(B, p,q),r[2],-l,9);} 
yBmp : : WriteBmpFile ( "optimized_coverage . bmp" , B) ; 
yKDT  ree : : Delete! ree ( N ) , delete [ ] B; 
for (int  i=0; i<10;++i)delete[ ]A[i] ; 

"^11  E  N  AUT  ^)GI^AI  L  •  LAS  T  PD  AT  E  D^l  S  D  U  L  2  0 
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OUTPUT: 


D=  4 

> 

D_bar=0. 262436 

y 

x=  0.183841 

y 

y=  0.031927 

1=  7 

y 

D_bar=0. 256637 

y 

x=  0.217818 

y 

y=  0.054442 

D=  8 

> 

D  bar=0. 241756 

y 

x=  0.324290 

y 

y=  0.017688 

D=12 

> 

D_bar=0. 239188 

y 

x=  0.303139 

y 

y=  0.112939 

D=13 

y 

D_bar=0. 238673 

y 

x=  0.321162 

y 

y=  0.061664 

D=14 

y 

D_bar=0. 209072 

y 

x=  0.597562 

y 

y=  Q.mniQ 

l_i 

II 

00 

y 

D_bar=0. 189872 

y 

x=  0.633593 

y 

y=  0.136296 

D=20 

y 

D  bar=0. 187058 

y 

x=  0.876571 

y 

y=  0.212578 

D=23 

y 

D_bar=0. 183453 

y 

x=  0.797085 

y 

y=  0.401234 

3=26 

y 

D_bar=0. 179227 

y 

x=  0.711760 

y 

y=  0.272756 

3=28 

y 

D_bar=0. 179046 

y 

x=  0.769481 

y 

y=  0.303093 

3=62 

y 

D_bar=0. 178932 

y 

x=  0.752994 

y 

y=  0.268364 

Fig.  8  Voronoi  diagrams  of  a  table  without  and  with  a  point  added  that  minimizes  D 


10.  Code  Summary 


A  summary  sheet  is  provided  at  the  end  of  this  report.  It  presents  the  yKDTree  namespace,  which 
contains  the  ColumnSort(),  NewTree(),  DeleteTree(),  and  NNSearch()  functions  and  the  NODE 
struct. 
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yKDTree  Summary 


y_kd_tree.h 


See  Yagerj  R.3.  "Working  with  KD-Trees 
Using  C++"  (ARL-TN-XXX) 


#ifndef  Y_KO_TREE_GUARD// 

#define  Y_KD_TREE_GUARD// 

namespace  yKDTree{//(^ _ 

template<class  T>void  CoiumnSoi^(//<=SORT  A  TABLE  BY  COLUMNS  (ASCENDING  ORDER) 

T**a,T**b,//< - - POINTERS  TO  STARTING  &  ENDING  ROWS  (FIRST  INDEX) 

T**t,//< . TEMPORARY  STORAGE  (SIZE  >=  b-a) 

unsigned  c){//< . THE  COLUMN  (2ND  INDEX)  TO  SORT  ON 

if (b-a<2)return; 

T**m=a+(b-a)/2]/*->*/ColumnSort(a,mjt,c),ColumnSort(mjb,tjC); 

for(T**u=fn, **l=a, j<b-a+tj )*j++=u<b&&(l>=m| |(*l)[c]>(*u) [c] ) ?*u++ :*!++; 

while (a<b)*a++=*t++; 

}//~^^~YAGENAUT@GMAIL.C<»1 - '-~LAST'-UPDATED~15:UL2014 - 

template<class  T>struct  NODE{//<================«==============A  KD-TREE  NODE 

T*rj//< . . . POINTER  TO  A  ROW  IN  A  TABLE 

NODE*a,*bj//< . POINTERS  TO  SUBNODES 

int  k;//< . NODE  LAYER 

};//~^~YAGENAUT@GMAIL.C0M~'-'-'— - - ~~~~~~^~~~~^lAST~UPDATED~153UL2014~~~— 

template<class  T>NODE<T>*NewTnee(//<==========================CREATE  A  KD  TREE 

T**a,T**b,//< . . . POINTERS  TO  STARTING  AND  ENDING  ROWS 

int  K,//< . NUMBER  OF  INDEPENDENT  DIMENSIONS 

int  k=0){//< . NODE  LAYER  (SET  BY  RECURSIVE  ALGORITHM) 

T**t=new  T* [b-a] ; /*->*/ColumnSort(aj  b,tj  k) j /*&*/delete[]tj 
T**m=a+(b-a-l)/2; 

NODE<T>*N=new  NODE<T>;/*<-*/N->r=*m,N->a=a-m?NewTree(a,m,K, (k+l)%K) :0, 
N->b=b-m-l?NewTree(m+lj  b,K, (k+l)%K) :0jN->k=kj 
return  N; 

}// - YAGENAUT@GMAIL.COM - — - — - LAST'-UPDATED~15:UL2014 - 

template<class  T>void  DeleteTree(//<==========================DELETE  A  KD  TREE 

NODE<T>*N){//< . . . THE  ROOT  NODE  OF  A  KD  TREE 

if( !N) return; /*else*/DeleteTree(N->a)jDeleteTree(N->b), delete  N; 

}// - YAGENAUT@GMAIL.COM - - - - - - - LAST'-UPDATED~15:UL2014 - 

template<class  T, class  U>void  NNSearch(//<=============NEAREST-NEIGHBOR  SEARCH 

const  NODE<T>*N,//< . . . THE  ROOT  NODE  OF  A  KD  TREE 

const  T*X,//< . POINTER  TO  SEARCH  COORDINATES 

T*&r,//< . POINTER  TO  NEAREST-NEIGHBOR  ROW  (CALCULATED) 

U&S,//< . SQUARED  DISTANCE  BETWEEN  COORDINATES  AT  X  AND  r  (CALCULATED) 

int  K){//< . NUMBER  OF  INDEPENDENT  DIMENSIONS 

NODE<T>*aj  *b; /*<-*/X[N->k]<N->r[N->k] ?a=N->a,b=N->b: (a=N->b,b=N->a); 
if (a)NNSearch(ajX,r jSj  K); 

U  s=0;/*<-*/for(int  k=0;k<K;++k)s+=(X[k] -N->r[k] )*(X[k] -N->r[k] ) ; 
if(s<S)r=N->rjS=s; 

if((s=X[N->k]-N->r[N->k])*s<S&&b)NNSearch(b,XjrjS,K); 
}//~^^~YAGENAUT@GMAIL.C0M~~'- - - - - - - - '-LAST~UPOATED~153UL2014— ~ 


#endif 


ColumnSortO  Performance 


#include  <algorithm>// . stable_sort  ( ) 

#include  <cstdio>// . printf() 

#in elude  <ctiffle>// . clock(),CLOCKS_PER_SEC 

#in elude  "y_kd_tree. h"// . yKDTree 

#include  "y_random.h"// . yRandom 

inline  bool  Compare(//<======HELPER  FUNCTION  FOR  THE  std: :stable_sort()  FUNCTION 

double*a,double*b){//< . POINTERS  TO  COMPARISON  ROWS 

return  a[0]<b[0];// . note  column  number  is  fixed  at  0 

- YAGENAUT^MAIL.COM—— - - - - - - ~LAST'-UPDATED~15:UL2014 - 

int  main(){//<==============MEASURE  THE  PERFORMANCE  OF  THE  ColumnSort()  FUNCTION 

const  int  N=l<<14jM=1000000;//. . . .max  #  of  rows,  number  of  iterations  per  test 

unsigned  I[625];/*<-*/yRandom: :Initialize(Ijl);//. .. .state  of  Mersenne  twister 

double  Sjt,*T[N],*A[N];/*<-*/for(int  i=0;i<N;++i)A[i]=new  double[l]; 


printf ( " 
row 
count 


I  std: :stable_sort()  | 

-I- 


time 

(s) 


I  Z[m/2][0]  I 
I  avg.  I 

-I . I- 


time 

(s) 


ColumnSort()\n" 

. \n" 

I  Z[m/2][0]\n" 

I  avg.\n" 

.| - \n");//table  header 


for(int  m=2;m<=N;m*=2){ 

s=0jt=clock() ,yRandom: : Initialize(1, 1) ; // . begin  stable_sort( )  test 

for(int  k=0;k<M;++k){ 

for(int  i=0;i<m;++i)A[i][0]=yRandom: :RandU(I,0,l); 
std : : stable_sort(A,A+m, Compare) , s+=A[m/2] [0]; } 
printf("%7d  |%9.3f  |%10.7f  | " ,m, (clock( )-t)/CLOCKS_PER_SEC,s/M); 

s=0jt=clock()  ,yRandom: :  Initialize(1, 1) ;  /  / . begin  ColumnSortO  test 

for(int  k=0;k<M;++k){ 

for(int  i=0;i<m;++i)A[i][0]=yRandom: :RandU(I,0,l); 
yKDTree: :ColumnSort(A,A+mjT,0), s+=A[m/2] [0] ; } 
printf("%9.3f  |%10.7f\n", (clock()-t)/CLOCKS_PER_SEC,s/M);} 
for (int  i=0; i<N;++i)delete[ ]A[i]; 

}// - YAGENAUT@GMAIL.COM~~'- - — - LAST~UPDATED~15:UL2014~~— - 


NNSearch ( )  Performance 


#include  <cstdio>// . . 

#include  <ctime>// . . 

#include  "y_kd_tree.h"// . . 

#include  "y_random.h"// . . 

templatecclass  T>T  NNSearchExhaustive(//<: 


. printf ( ) 

. clock(),CLOCKS_PER_SEC 

. yKDTree 

. yRandom 

^EXHAUSTIVE  NEAREST-NEIGHBOR  SEARCH 


J**a,J**b,//< . . . POINTERS  TO  STARTING  AND  ENDING  ROWS 

T*X,//< . POINTER  TO  SEARCH  COORDINATES 

T*&r,//< . POINTER  TO  NEAREST-NEIGHBOR  ROW  (CALCULATED) 

int  K){//< . NUMBER  OF  INDEPENDENT  DIMENSIONS 


double  S=0;/*<-*/for(int  i=0; i<K;++i)S+=(a[0] [i]-X[i])*(a[0] [i] -X[i] ); 
for(r=*a; a<bj++a){ 

double  s=0; /*<-*/for(int  i=0; i<K;++i)s+=( (*a) [i] -X[i])*((*a)[i]-X[i]); 
if(s<S)r=*a,S=s;} 
return  S; 

}// - YAGENAUT@GMAIL.COM—'— - - - LAST'-UPDATED'-15:UL2014 - - 

int  main(){//<================MEASURE  THE  PERFORMANCE  OF  THE  NNSearch()  FUNCTION 

const  int  N=1«14,M=10000000,K=2; 

unsigned  I[625];/*<-*/yRandom: :Initialize(Ijl);//. .. .state  of  Mersenne  twister 
double*A[N] jB[N][K];/*<-*/for(int  i=0;i<N;++i)A[i]=B[i]; 
for (int  i=0; i<N;++i)fon(int  k=0;k<K;++k)A[i] [k]=yRandom: : RandU(I,0,l) ; 
printf("  I  NNSearchExhaustive() I  NN2DInterp()\n" 

"  row  I - I - \n" 

"  count  I  time  |  x  |  time  |  x\n" 

"  I  (s)  I  avg.  I  (s)  I  avg.\n" 

” - I - I - I - I - \n" );//.  .table  header 

for(int  m=2;m<=N;m*=2){ 

double  s=0,t=clock();/*&*/yRandom: :Initialize(I,l);//. . . NNSearchExhaustive( ) 
for(int  k=0;k<M;++k){ 

double  X[K] ; /*<-*/for(int  i=0; i<K;++i)X[i]=y Random : : RandU(I,0, 1) ; 
double*r; /*<-*/NNSear ch Exhaustive (A, A+mjX, r, K) ; 
s+=*r; } 

printf ("%7d  |%8.3f  |%9.6f  | " ,m, (clock( ) -t)/CLOCKS_PER_SEC,s/M); 

s=0, t=clock(), yRandom:  :Initialize(I,l);// . begin  NNSearch()  test 

yKDTree : : NODE<double>*R=yKDTree : : NewTree(A,A+m,K); 
for(int  k=0;k<M;++k){ 

double  X[K] ; /*<-*/for(int  i=0; i<K;++i)X[i]=y Random : : RandU(I,0, 1) ; 
double  S=1E9; 

double*r;/*<-*/yKDTree: :NNSearch(RjX,r,SjK); 
s+=*r; } 


yKDTree : : DeleteTree(R) ; 

printf ( "%8. 3f  |%9. 6f \n", (clock() -t)/CLOCKS_PER_SEC, s/M) ; } 
}//...w~..^~YAGENAUT@GMAIL.COM'-'-'-~~'— - - - - - '-'-LAST'-UPDATEO'-15:UL2014'-'-— 


K=2 


K=4 


K=6 


Example  -  Creating  a  Simple  Voronoi  Diagram 

#include  <cmath>// . fabs() 

#in elude  "y_bmp.h"// . yBmp,<cstring>{memcpy()} 

#in elude  "y_kd_tree .  h "// . yKDT ree 

inline  void  Rainbow(//<========================================RAINBOW  COLOR  MAP 

unsigned  char  C[3],//< . OUTPUT  COLOR  (CALCULATED) 

double  x,//< . VALUE  FOR  WHICH  A  COLOR  WILL  BE  CALCULATED 

double  min, double  max){//< . MINIMUM  AND  MAXIMUM  SCALED  VALUES 

if(x<min){C[0]=C[l]=C[2]=0;/*&*/return;}// . set  too  small  values  to  black 

if(x>max){C[0]=C[l]=C[2]=255;/*&*/return;}// . set  too  large  values  to  white 

x=(l-(x-min)/(max-min))*8j// . remap  x  to  a  range  of  8  to  0 

C[0]=int((3<x&&x<5| |x>7  ?-fabs(x/2-3)+l. 5: 5<=x&&x<=7?l:0)*255) ; // . blue 

C[l]=int((l<x&&x<3i i5<x&&x<7?-fabs(x/2-2)+1.5:3<=x&&x<=5?l:0)*255);// _ green 

C[2]=int((  x<li i3<xS&x<5?-fabs(x/2-l)+1.5:l<=x&&x<=3?l:0)*255);// . red 

}// - YAGENAUT@GMAIL.COM - - - — - LAST'-UPDATED'-15:UL2014 - — 

int  main(){//<=========================CREATE  A  VORONOI  DIAGRAM  USING  NNSearch() 

double*A[5],B[5][3]={4,l,0  ,  4,3,1  ,  6,2,2  ,  2,4,3  ,  8,4,4}; 
for(int  i=0;i<5;++i)A[i]=B[i]; 

yKDTree: :NODE<double>*N=yKDTree: :NewTree(A,A+5,2); 

int  n=1000;// . image  size  will  be  2n  x  n  pixels 

unsigned  char*I=y8mp: :NewImage(2*n,n,255),BLACK[3]={0}; 
for(int  p=0;p<2*n;++p)for(int  q=0;q<n;++q){ 
double  X[2]={p*10./(2*n-l),q*5./(n-l)>; 
double  S,*r;/*<-*/yKDTree: :NNSearch(N,X,r,S=lE9,2); 
if(S<.002)memcpy(yBmp: :GetPixel(I,p,q)jBLACK,3); 
else  Rainbow(yBmp: :GetPixel(I,p,q), r [2] ,  -1. 5,5 . 5); } 
yBmp: :WriteBmpFile(" voronoi. bmp", I); 
yKDT  ree : : Delete! ree(N) , delete [ ] I; 

}//...w....^~YAGENAUT@GMAIL.COM'-'-'-— - - - - - LAST'-UPDATEO'-15:UL2014'-'- - 
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